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Abstract Species do not merely evolve, they also coevolve with other organisms. 
Coevolution is a major force driving interacting species to continuously evolve ex- 
ploring their fitness landscapes. Coevolution involves the coupling of species fit- 
ness landscapes, linking species genetic changes with their inter-specific ecological 
interactions. Here we first introduce the Red Queen hypothesis of evolution com- 
menting on some theoretical aspects and empirical evidences. As an introduction 
to the fitness landscape concept, we review key issues on evolution on simple and 
rugged fitness landscapes. Then we present key modeling examples of coevolution 
on different fitness landscapes at different scales, from RNA viruses to complex 
ecosystems and macroevolution. 



1 Introduction: the Red Queen 

Coevolution pervades evolutionary change on multiple scales. It is not exaggerated 
to say, changing a little the classical Dobzhansky's statement, that nothing makes 
sense in biology except in the light of coevolution. Darwin himself recognized this 
when referring to what he called the entangled bank (9): "It is interesting to contem- 
plate an entangled bank, clothed with many plants of many kinds, with birds singing 
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on the bushes, with various insects flitting about, and with worms crawling through 
the damp earth". Indeed, ecosystems need to be seen as collectives of interacting 
species whose evolutionary fate is necessarily intermingled in complex ways. Such 
complex networks pervade ecosystems and their evolutionary dynamics [48 1. 

The rationale of the previous statements is simple. As any species changes over 
time, it inevitably triggers co-evolutionary responses in those partners directly af- 
fected by their interdependencies. A prey running faster than its predator will require 
changes in the later to cope with the change. Running fast is one option, hiding in 
the appropriate place another. If we move ourselves into the microcosmos of host- 
pathogen interactions, including cells and viruses (or bacteria) as main examples of 
coevolving players, a faster rate of change in the parasite might need a faster re- 
sponse from the host. Changing attributes might be an unavoidable consequence of 
entangled ecosystems and not changing might possibly be lethal. This view matches 
a dynamically unstable scenario where changes keep happening all the time and, as 
in Lewis Carroll's Through the looking glass, species need -as Alice does- to con- 
stantly run to remain in the same place. Such picture was early supported by Leigh 
Van Valen's work, and is known as the Red Queen hypothesis. 

The Red Queen model was introduced by Leigh Van Valen in 1973 |87|. It was 
conceived as a theoretical explanation for the observation that the extinction prob- 
ability of a species is approximately independent of its length of existence Il8~7l l4ll. 
Accordingly with this view, Van Valen observed that the vast majority of taxonomic 
groups analyzed displayed exponentially decaying survivorship curves. This result 
implied constancy in the probability of extinction of the taxa, regardless of their 
previous duration. That is, both data from the fossil record and from extant species 
suggested that a given species may disappear at any time, irrespective of how long 
has already existed. This unexpected phenomenon, termed the Law of Constant Ex- 
tinction, can be formulated, in a simple way, as follows. If N(t) indicates the number 
of species at a given time and we follow their presence over time (ignoring other 
events) we would observe an exponential decay law, namely: 



where 8(t) indicates a time-dependent extinction rate. If Nq is the original number, 
this differential equation is easily solved, and gives: 



Despite the seemingly obvious assumption that 8 depends on t, the surprising ob- 
servation is that the observed curves fit very well a constant decay rate 8, i.e., a 
solution: 



where 8 is the extinction probability of a species (per millions of years, Myr). This 
law is essentially correct on average, despite the fine-scale pattern is much more 
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Fig. 1 Extinction and the fossil record of life, (a) Successive cohort survivorship curves for 2,316 
extinct marine families during the Phanerozoic (redrawn after |57|). Notice that, together with 
an average exponential tendency to decline (which would give straight lines in a linear-log plot), 
discrete and punctuated events are found in the dynamics of survivorship. This pattern is well 
appreciated in (b) where we display the number of ammonoid genera surviving over geologic time. 



episodic ll56l . as depicted in figure [TJ a). Here we display the surviving sets of fam- 
ilies found in the marine fossil record (the so called seudocohorts) through time. A 
roughly exponential decay can be identified, together with sharp extinctions events 
(see also figure[TJb)). 

Our intuition, guided by Darwin's theory of natural selection, would have ex- 
pected species within any group to become longer lived along time: if adaptation 
improves species progressively through time, a decreasing probability of extinction 
should be expected. That is, older species should last longer. However, careful anal- 
ysis revealed that species of modern mammals are likely to become extinct as were 
their ancestors living 200 Myr ago 0). If evolution leads to improvement through 
adaptation, why modern mammals have the same extinction probabilities as their 
ancestors? Van Valen's interpretation is simple but counterintuitive: species do not 
evolve to become any better at avoiding extinction. Van Valen suggested that con- 
stant extinction probability would arise in an always changing biotic community, 
with species continually adapting to each other's changes. The name for his conjec- 
ture alludes to the Red Queen's remark in Lewis Carroll's Alice Through the Looking 
Glass: "here, you see, it takes all the running you can do, to keep in the same place". 
Van Valen's view of evolution is that species change just to remain in the evolution- 
ary game and extinction occurs when no further changes are possible. Actually, the 
Red Queen hypothesis is profoundly Darwinian, in that it puts emphasis mostly on 
biotic interactions rather than on abiotic factors [29|. Van Valen further elaborated 
this concept in much more detail in subsequent articles, showing that his theory was 
compatible with the classic population-genetic view of species evolution Il85ll86l . 

To test the plausibility of the Red Queen hypothesis, Maynard Smith and co- 
workers (see ||79l and references therein) developed a theoretical model describing 
continuous (co-)evolution of species in a constant environment. Such model consid- 
ered a fixed number of S interacting species, defining some fitness measure 0, and 
a maximum fitness 0* was supposed to exist for each species in a given fixed, ex- 
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ternal biotic environment. At a given time, the fitness 0, and the maximum 0* took 
different values, and each species "tried" to reduce the so called lag load, defined 
as: 

<t>i - </>* 
L ; = ^-^; i=l,...,S. 

0i 

If Pij is the change in the lag load L, due to a change in Lj, then a mean-field equa- 
tion for the average lag load (L) = can be derived. This is done by first 
separating, for each species, changes due to "microevolution of coexisting species" 
from those linked with its own microevolution [79]. The whole equation for the 
lag load variation in a given species is: 8L{ = 5 C L, — 8„Li, which simply says that 
it typically increases due to changes in the other species and decreases due to mi- 
croevolutionary changes in the species under consideration i.e., 5 C L, is the increase 
in the lag of the /th species caused by evolutionary changes in others, and 8 g Lt is 
the reduction in lag caused by changes in species i itself. This can be written in the 
following way, 

s 

8 Li = £ Pi jSgLj - 8 s Lj, 
;"=i 

where j3y (with /?,-; = 0) is the increase in L, due to a (unit) change in Lj, Assuming 
that most species are close to their adaptive peaks, any evolutionary change in one 
species will have a deleterious effect on the other species. The deterministic, time- 
continuous equivalent model can be formulated with: 

dL s 

ai j=l 

By taking the average in both sides of Eq. ([TJ, the evolution of the average lag 
load is given by: 

Assuming now that ki = k for all i = 1, ...,S, the average lag load equation can be 
written as: 

7=1 

and it has a steady state solution if f, = 1 for all j = In other words, if: 

r = £j3 l7 = l; Vj. 

Otherwise, it can be shown that (L) will decrease (increase) for F < 1 (T > 1). The 
previous identity is telling us that the equilibrium state of this system is reached 
through a balance between the reduction of the individual lag load of each species 
and the increases due to coevolutionary changes in the remaining partners. The main 
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result of this model is that evolution of species proceeds at an approximately steady 
rate even in the absence of external or environmental changes [79|. At the end of this 
chapter we will present a dynamic model of Red Queen dynamics where evolving 
networks interactions are made explicit. 

As we previously discussed, the Red Queen hypothesis provides a plausible ex- 
planation of the fossil data record, but it turned to have more implications. For in- 
stance, one suggested implication of the Red Queen hypothesis is that coevolving 
pathogens may facilitate the persistence of outcrossing despite its costs. Specifi- 
cally, coevolutionary interactions between hosts and pathogens might generate ever- 
changing conditions and thus favor the long-term maintenance of outcrossing rel- 
ative to self-fertilization [2] or asexual reproduction l35l l26l (see also |27| for a 
review). Outcrossing (mating between different individuals) involves the introduc- 
tion of unrelated genetic material into a breeding line, thus increasing genetic diver- 
sity. The previous statements are supported by evidences from nature. For instance, 
experimental studies on the coevolution of a nematode with a bacterial pathogen 
[49 1 revealed that the action of parasites caused an increase of outcrossing in mixed 
mating populations. Interestingly, these experiments also revealed that coevolution 
with the pathogen caused extinction in populations without outcrossing, whereas 
outcrossing populations persisted through reciprocal coevolution. Studies in natu- 
ral snails populations also revealed that sexual reproduction is more common when 
parasites are abundant and adapted to infect local host populations l43l l40l . Co- 
evolution and Red Queen dynamics were also identified for the crustacean genus 
Daphnia and its parasites in pond sediments ifTUl . 



2 Red Queen on a Lattice: a toy model 

Before we get into the more formal approaches taken to describe and simulate the 
evolution and coevolution of species on fitness landscapes, let us consider a simple 
toy model that illustrates the basic idea behind van Valen's metaphor. Imagine a 
world where our species can move on the surface of a sphere. To makes things 
simpler, consider a discretized surface, like a mesh covering the spher^j] 

To simulate such a system we used the so-called cellular automata (CA) models 
(32). CA models are a common tool to investigate interacting agents in a physical 
space, which, for our case, will be a surface. Each point in this mesh is a site, which 
can be either empty or instead occupied by an individual of a given type. Let us start 
with a simple "ecosystem" formed by a species exhibiting two phenotypic traits. Let 
us indicate as E — {0, 1} the two possible "genotypes" which can be understood as 
two alleles of a given gene. 

In our idealized model, and 1 are the only two genotypes, each one associ- 
ated to a set of parameters defining the underlying phenotype. For simplicity, we 

1 Specifically, we start from a lattice, whose surface has been discretized using a mesh, and then 
we perform a projection of this mesh on a surface by properly deforming the initial coordinates 
using a so called icosahedron-based pixelization. For details, see 1 80 1 . 
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Fig. 2 Spatial competition dynamics in a single-species model with two mutants. Starting from 
a small set of occupied sites with and 1 genotypes (white and black, respectively) scattered on 
an empty (gray) landscape, we display spatial patterns for different simulation times, T. After a 
few steps, local populations start to grow, but competition is still weak. Once the two populations 
occupy enough space, competition starts to be effective but no global exclusion occurs. Instead, the 
two populations coexist by expanding locally over spatial domains that appear homogeneous. After 
a long time, the system is rather stable. Although structures keep changing their boundaries, the 
global picture remains the same, with a spatial landscape displaying large homogeneous patches. 
Here we used [1 = 10~ 4 , 8h = 10~ 3 and r H = 0.35. 

will consider completely symmetric sets. In other words, we have a neutral change 
when moving from to 1 and viceversa. These transitions occur proportionally to 
mutation probability ji £ [0,1]. 

We can use the previous system to simulate host-parasite interactions with match- 
ing alleles (MA, see Section 4 for further details) interactions. Our CA is thus given 
by a two-dimensional state space, Q(iJ), with spatial coordinates The states, 
S, of the automaton at time t are given by S(i,j;t) G E = {Hk,Pk,E} , where H and P 
denote, respectively, hosts and parasites defined as 1-bit strings (i.e., with k = 0, 1). 
E indicates empty sites in the state space Q. 

CA models include dynamics by means of the state-transition rules, which deter- 
mine the changes of the states according to the current states in a given site together 
with its neighboring states. For this particular system we used the so-called Moore 
neighborhood, which considers the eight nearest neighbors. The model rules are 
summarized in Box 1. 
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Fig. 3 Red Queen dynamics in host-parasite (prey-predator) coevolution. Here we have added to 
the only-host system displayed in the previous figure an evolving 1-bit parasite. Parasites propagate 
through space provided they find the same host genotype. Here two levels need to be considered: 
the spatial framework (a) defined by the local arrangement of individuals and the required geno- 
typic matching, (b) Schematic diagram of the two populations (upper and lower layers) with back 
and forth mutations (arrows) between genotypes and the requirement of allele matching (vertical 
lines). The plots in (c) show the host populations (as before) which now keep changing. Simi- 
lar plots would be observed (although having less dense patches) for parasites. Here we used, as 
before, n„ = 1(T 4 , S H = 1(T 3 , r H = 0.35; and fi P = 1(T 4 , 8 P = 0.2, and r P = 0.1. 



Box 1. State transition rules of the CA model for host-parasite matching allele 
interactions with 1-bit genotypes. 

1 . Death: hosts and parasites decay with probabilities 8h and 8p respectively, 
according to (recall k = 0, 1): 

H k %E, 
Pk^E. 

2. Birth of new hosts: hosts can replicate with probability r#, without and 
with mutations, following the reactions: 

H k ^H k +H^ kl 
where jj.^ is hosts mutation rate. 
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3. Predation: Predator genotypes predate on hosts same genotype (assuming 
a perfect MA interaction), with rate rp. In other words, they can reproduce 
only under the presence of the same host genotype in the neighborhood. If 
such condition is met, the new reactions follow: 

n + tlk > Lrk , 

Pk+H k ^P k +P^ kl 

introducing again mutational changes proportionally to parasites mutation 
rate jlp. Because of the H-P interaction, hosts experience a parasite-driven 
mortality. 

What is the dynamics of this simple model when no parasites are present? In 
figure [2] we see an example of the time evolution of this model, starting from an 
initial condition where we scatter a small population of each genotype over the 
sphere. The empty, available space is indicated in gray, whereas the two alternative 
genotypes are shown as black and white squares. After a short transient, where both 
variants expand with no special constrains, available sites become scarce and the 
expanding patches grow and develop rugged boundaries. Such pattern stabilizes in 
the long run, where we observe large domains of each class. This phenomenon is due 
to the local exclusion of our identical competitors [7 1 1 that allows global coexistence 
to occur. 

The previous scenario shows that competing populations end up in a predictable 
community structure with no further (global) changes, However, when an additional 
component -parasites (or predators)- is added in the same system, it immediately 
triggers the emergence of an unstable dynamical state. This is shown in another 
example in figure [3] where we again represent the host populations, now starting 
with the same condition as in figure [2] but adding also some randomly distributed 
predators. 

At any step, the spatial distribution changes rapidly and complex waves of ex- 
pansion and contraction, affecting both genotypes, are observable. Sometimes, the 
extinction of the parasites returns our system to the previous conditions without par- 
asite (and a spatial segregation pattern). This occurs for example when the mutation 
rate of the parasite is too small or its death rate too high. Sometimes, the pressure of 
the parasites is so strong that they cause the extinction of the hosts and the eventual 
collapse of the whole host-parasite system. The interesting dynamics occurs when 
parasites are able to reliably match their preys and reproduce at a reasonable pace. 
Similarly, mutation rates need to be high enough to react to depleted host popu- 
lations and at the same time allow for a conservation of genetic information, thus 
avoiding undesirable drift (see Section 5.1). In other words, intermediate rates of 
parasite pressure end up in Red Queen evolution. 

Here the system constantly changes as a consequence of the hide-and-run effect 
induced by the parasitic species. Each time a parasite finds a suitable host sharing 
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the same bit, it replicates and will keep expanding provided that it finds additional 
hosts in the neighborhood. The end result is a spatial pattern of propagating patches 
and constantly changing distributions of the two available genotypes. 



3 Fitness Landscapes 

Adaptive or fitness landscapes are a very useful tool and an important concept in 
evolutionary biology. They are used to map, represent or visualize the relation- 
ship between genotypes (or phenotypes) and reproductive success or fitness. Fit- 
ness landscapes were first introduced by Sewall Wright and were then extended by 
other authors (see [38, 54] and references therein). Fitness landscapes assume that 
each genotype has a well-defined fitness value, which is represented with a height 
or peak in the landscape (see figures [5] and [8]). A landscape simply means a single- 
valued scalar function, F(x), of the state or configuration x of a system. The variable 
x typically has very many dimensions, and thus may be often written like a multidi- 
mensional vector, as a set of N components xf. 

x= (xi,x 2 ,...,xn). 

The term landscape is inspired from the geographic landscapes in which the height 
h above sea level is a simple function h = F(x,y), of the two-dimensional location 
x = (x,y). In the field of biology, fitness landscapes are generically representing the 
fitness of a given biological entity as a function of its genotype or phenotype. As 
biological entity we referee to a given organism or a to particular macromolecule or 
cell of that individual. 

Fitness is a relative measure, since it may depend on the environment and on 
other interacting organisms ll54l . Fitness can be given by several properties, or by a 
combination of them. For instance, we can use replication or reproduction success 
as a measure of fitness. Properties like infectivity, migration capacity, ability to co- 
operate, among others, can also define a fitness which may facilitate survival and 
adaptation. As Jacob ll33l[34l stated, adaptation typically progresses through small 
changes involving a local search in the space of possibilities (e.g., sequences space). 
The paradigm is a hill-climbing process via fitter mutants which "move" towards a 
global or local optimum in the fitness landscape (see figure|4|A)). The hill-climbing 
framework was originally proposed by Wright ||90ll9T1 l, who introduced the concept 
of space of possible genotypes. Under this framework, each genotype has a given 
fitness, being the distribution of fitness values over the space of genotypes a fitness 
landscape. Depending upon the distribution of fitness values, the fitness landscape 
will become more or less mountainous. The behavior of an adapting population will 
depend on how rugged the fitness landscape is, on the size of the population, and 
on the mutation rate which moves a population from one genotype to another in se- 
quence space. The motion of a population over a fitness landscape also depends on 
whether the population reproduces asexually or sexually. The latter reproduction in- 
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volves mixing of genotypes which can involve to reach more distant points in fitness 
landscape, compared to asexual reproduction |38|. 

The fitness landscape concept has been widely used in both evolutionary and 
coevolutionary biology. Organisms do not merely evolve, they coevolve with other 
organisms. As a difference from evolution, that can be roughly characterized as an 
adaptive search on a "potential surface", in coevolution there may typically be no 
such potential surface, being the process far more complex [38]. Actually, in coevo- 
lutionary processes the adaptive landscape of one organism can deform and heave 
as the other organisms make their own adaptive moves. Under this perspective, one 
can interpret coevolution as both dynamical and evolutionary processes occurring 
in coupled fitness landscapes (see figures [3jb),|4]|5][8]and[T0ja)). 

For the sake of clarity, before accounting for coevolution, we will introduce, 
in Sections 3.1. and 3.2., some information about evolution on fitness landscapes. 
First, we will describe evolutionary phenomena in simple fitness landscapes also 
presenting a theoretical body to model evolution in these types of landscapes. Then, 
we will extend our explanations to more complex, rugged fitness landscapes. From 
Section 3.3. onwards (together with the Introduction Section above) we will strictly 
focus on coevolutionary phenomena. As an example of coevolution at small scales, 
Section 3.4. includes the view of RNA virus evolution from the perspective of the 
Red Queen hypothesis. The remaining sections will deal with some examples and 
models in higher biological scales, from complex ecosystems to macroevolution. 



3.1 Simple versus coupled fitness landscapes 

Models on evolution have considered different theoretical and computational frame- 
works to characterize several levels of complexity. One of the most successful ap- 
proaches to address evolutionary phenomena (as well as coevolution as we will 
discuss later) on fitness landscapes is given by the digital genomes approach ETI . 
Under this approach, using as an example the evolution of RNA genomes, we can 
develop a mapping between RNA sequences (defined as a chain of nucleotides in- 
volving a four-letter alphabet £2) and binary sequences, according to: 

J? : Q = {U,G,A,C} — > E = {0, 1}. 

Alternatively, one can use another Boolean representation using spins instead of 
bits: 

: Q = {U,G,A,C} — > E = {+1,-1}. 

Both approaches are equivalent because the mapping has the same nature. How- 
ever, the spins approach exploits some advantages of considering "up" and "down" 
configurations to describe the microscopic dynamics. 

Let us define the ith string of the population, S; = (Sn , .. . , Sj V ), as digital genomes 
(i.e., sequences of purines and pyrimidines only incorporating the linear information 
encoded by the string) of length v. In order to determine the functional relevance of 
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these sequences, we need to map them to a sequence-fitness measure, which can be 
defined as: 

/:S,-: L v ./(S,i. 

This functional relation can be generically divided into two types: (i) the different 
bits of a string play an independent role in fitness; or (ii) some bits of the string influ- 
ence others in a nontrivial way. Case (ii) corresponds to what is known as epistasis. 
Epistasis occurs when the phenotypic effect of a mutation depends on the presence 
of other mutations in the genome. Epistasis becomes especially important in highly- 
compacted genomes that are expected to contain multifunctional proteins or over- 
lapping genes. In this sense, epistasis has been studied and characterized for RNA 
viruses, both experimentally and theoretically (see [ 18 1 for a review and |64 65|)- In 
a more general way, epistatic interactions play an important role in evolutionary ge- 
netic systems almost whenever multi locus genetics matters and plays a central role 
in the evolution of genetic systems such as sex and recombination, ploidy, genomic 
segmentation and modularity, genetic incompatibility and speciation, mechanisms 
of mutational robustness, mutational load for deleterious mutations through genetic 
drift, and the rate of adaptive evolution 1 1 2 1 . 

The digital genomes approach allows us to use an abstract, multidimensional 
representation of the potential set of states accessible to a v-bits digital genome. 
This set or space is given by a sequence space in the form of a hypercube, v = E v , 
which can provide, at low dimensions, some intuitions about the behavior of strings 
under selection-mutation pressures (see figure |4}. If only small mutation rates are 
considered, transitions between sequences will take place only involving nearest 
neighbors in sequences space, thus differing only in one bit. 

In general, for a given mutation rate, ji, two sequences S and S' will be obtained 
from each other with a given probability, given by: 

(S -> S') = ^«( s ' s ') (1 - jU )V-^(S,S') ) 

where c///(S,S') is the Hamming distance among the two sequences (i. e., the 
number of different bits), with: 

d*(S,S') = £(l-%,A 

i=l V ' 

where 5/j is Kronecker's delta with 8i j = 1 if i = j and = if i ^ j. Here 
can be interpreted in probabilistic terms: it is the probability of having exactly cLh 
differences between the two digital genomes. This function allows to introduce the 
dynamics associated to mutations as transition probabilities. For the spin mapping, 
the transition probabilities can be expressed as: 

^(S^S')=^exp(-j3£ SiSM, 
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Fig. 4 (A) Schematic representation of a fitness landscape with three peaks. Depending on the 
initial condition in the genotype space, the population will evolve towards different maxima due to 
mutations which result in motions in the landscape. (B) Fitness landscapes can also be represented 
in the sequences spaces (here using digital genomes). Such systems allow defining trajectories 
followed through string evolution. Four standard cases of 3-dimensional sequences space (where 
the size of the nodes denotes each string fitness value) are also displayed: (a) flat, (b) Fujiyama, (c) 
Swetina-Schuster (single peak), and (d) rugged fitness landscapes. 



where the /3 term, defined as /3 = log(/i/(l — ju))/2, can be interpreted in terms of 
a temperature, being jV a normalization constant. 

Once we define the fitness function associated to each vertex of the hypercube, 
we can characterize the dynamics. 

If N(S,t) indicates the fraction of strings having a given sequence S £ E v at time 
f, Eigen's formulation ifTTIl describes the population dynamics as a set of nonlinear 
differential equations, given by: 

^^=£W M (S'^S)/(SW ( 2 ) 

The first term on the right-hand site of Eq. (|2]i corresponds to positive contri- 
butions to the abundance of S due to mutation transitions from other strings of the 
sequence space. The second term includes all the reverse events leaving the node 
occupied by S. Figure^B) illustrates the information described in Eq. Q for 3-bits 
strings (i.e., v = 3). The nodes of the hypercube indicate the population size and the 
three potential transitions from 000 to other strings differing in one, two or all bits 
are indicated by arrows of different colors (see [ 18 1 for the extension of the previous 
results to discrete dynamical systems). 
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Simple fitness landscapes can be defined from our previous definitions. Roughly, 
a sequences space is a discrete space including all sequence combinations of a 
given sequence, which are connected to neighboring sequences differing in one bit. 
This space thus results in a set of nodes or vertices (sequences) which are con- 
nected by single-point mutations (see figure[4|a-d)). For a given sequence of length 
V, the dimension of the sequences space is given by Jf v . For DNA or RNA se- 
quences, the total number of nodes in the sequences space will be 4 V , which results 
in an astronomic number. As we saw before, this inherent complexity can be re- 
duced using digital genomes defined as bit strings or as "up" and "down" spins. 
The bit-strings approach allows one to simulate the processes of (co-)evolution 
and selection under different types of rules or interactions describing different pro- 
cesses in biological systems. For instance, in cancer dynamics |69 1, in RNA viruses 
llS4ll70ll4Tl[T8ll65ll66ll68l , and, in the context of coevolution in matching-alleles 
dynamics (see Section 4), among others. As a simple, illustrative example, figure|4] 
shows four different simple fitness landscapes. The sequences space in figure |4| a) 
is a flat fitness landscape, where all the sequences have the same fitness /o, ac- 
cording to J^(Si) = fo- If the sequences space has a least fit sequence (e.g., 0) 
and then the fitness increases at increasing number of mutations, we have the so- 
called Fujiyama fitness landscape, shown in figure |4|b), being its fitness given by 
J$?(Si) = fy + YX=\ $ik ( see [55 1 for the application of the Fujiyama landscape to 
RNA viral populations). Another widely studied case is the so-called single-peak 
fitness landscape, which is shown in figure |4jc). For this landscape, the fitness can 
be defines as JT(S/) = f 8s h i +/i (1 - 8s u i), with f > f x (see l68l l63l l46l for 
some examples and applications of the single-peak fitness landscape). 



3.2 Evolution on rugged fitness landscapes 

The sharp, single-peak fitness landscape cited at the end of the previous section 
defines an extreme in a hierarchy of models introducing different levels of depen- 
dencies among genes. A different approximation deals with landscapes in a much 
more general way, by allowing them to display a given number of local maxima 
generating a mountainous landscape. Together with the simple fitness landscapes 
we display a rugged fitness landscape in figure^d). For this landscape, where each 
sequence has a different fitness resulting in as many peaks as sequences, the fitness 
can be given by ^(S,) = hY%=i%ik> being ^ e fc[0,l], a random number. The 
best known model for the evolution on rugged fitness landscapes is Kauffman's NK 
model [38 , 37 1, which is also defined on a hypercube. It was originally proposed as 
a representation of haploid genomes involving two alleles per locus with additive 
contributions to fitness from different loci. 

The NK model is a simple model of random epistatic interactions. In this model 
N is the number of parts of a system (e.g., genes in a genotype or amino acids in 
a protein). Each part makes a fitness contribution which depends upon that part as 
well as upon K other parts among the N. Thus, K reflects how richly cross-coupled 
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the system is indicating how many other genes influence a given gene i.e., the rich- 
ness of epistatic interactions among the genes. Such a model, under parameters al- 
teration, generates a family of increasingly rugged multi peaked landscapes. Once 
again a fitness function is introduced, f = f(Sn,...,Si V ), and changes in the traits 
are assumed to occur by means of single, one-bit steps (i.e., single-point mutations). 
This single -chain events are consistent with our assumption of small mutation rates. 
In this way, a given string obtained by inaccurate replication allows to perform a 
random adaptive walk from a given node towards one of its v nearest neighbors 
if this leads to an increase in fitness. A direct consequence of this process is that 
once a local peak is reached, no further changes are allowed to occur. This is com- 
pletely different from the assumptions made above, which assume the presence of 
a preferred sequence around which other sequences have a lower fitness value. In 
the context of NK landscapes, a local peak is very simply defined: if all nearest 
neighbors in the hypercube are less-fit, we have a fitness local maximum. 

How can we construct a system displaying a NK landscape? Kauffman suggested 
a simple approach using fitness tables: for each element Sy, if it is influenced by K 
other elements, each element contributes in an additive way to the overall fitness. 
In other words, if we consider the two-locus model and assume that a given locus 
i constributes to the global fitness associated to S by an amount /j(S) G [0, 1], the 
global fitness is given by the average value: 

f(s)=lifi(s). 

v (=1 

As K grows, the ruggedness of the landscape increases, since more complex inter- 
actions are allowed to occur. 

An interesting feature of the NK model is that, because of its simplicity, it al- 
lows the prediction of some evolutionary dynamical patterns. As an example, let us 
consider that fitness values are random and uncorrelated, i. e., f(Sn,...,Si V ) = 
where | £ [0, 1] is a random number with uniform distribution. This random fitness 
landscape has many local fitness peaks. This number Ml is very large: 

Mi(v) = vTT' 

and thus our digital sequences can get trapped in a very large number of optima. 
To see this, let us consider the number of neighbors of a given node and compute 
the probability that this node is a local maximum. The chance that it is the fittest 
among its v neighbors and itself, given the random choice of values, is simply Pi = 
l/(v+ 1). Since there are 2 V possible strings, the fraction of those who are local 
maxima is Mi(v) = 2 V P\. An extension of this model can be easily introduced by 
means of the so-called Fujiyama landscape (see previous Section), where a fitness 
function is defined now as follows: 

f( Si )= l -(l-s) k , 
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with k = v — Yj=i Su, being z = Lj=i (1 — sV a normalization factor. The parameter 
s weights the steepness of the peak. 

Similarly, the presence of epistatic interactions can be introduced using the 
sequence-dependent fitness: 



m) = 



i 



</«(s,-,s„) 5 



where | > defines the degree and type of epistasis (for | < 1 we would have 
antagonistic interactions whereas S, > 1 defines synergistic epistasis (see [18] and 
references therein). 



3.3 Coevolution on rugged fitness landscapes 

The previous sections have been developed to provide the reader with a broad frame- 
work that will be, from now on, extended to the subject of coevolution. As we al- 
ready focused on evolution on rugged fitness landscapes, let us start with coevo- 
lution on these landscapes (in following sections we will analyze coevolutionary 
dynamical models in simpler fitness landscapes). In the context of rugged fitness 
landscapes, the NK model can be modified to analyze evolution between many in- 
teracting species, by means of the so-called NKC model l37l . This model introduces 
a new parameter, named C, which denotes the number of couplings between differ- 
ent species (also represented as strings). Now, the fitness of the NK model needs 
to be modified in order to introduce the coupling effects where each trait receives 
inputs from other C other traits from different species. These traits are randomly 
chosen between the S species of the ecosystem. The NKC model includes three 
main parameters describing: (a) the number of traits required to characterize a given 
species (AO, (b) the number of so-called epistatic interactions among genes in the 
same species (K), and (c) the number of interactions among traits of different species 
(C), which introduce coevolution. 

Figure [5] illustrates this approach for N = 3 traits of two interacting species. In 
the figure, the local peaks are indicated by black circles. Each species is defined by 
a set of traits, which are coded by bit-strings. Such traits are connected among the 
different species. In figure |5j a), species 1 is not located in a local peak. As a result 
of an adaptive walk, it will reach the local peak by mutation. However, as a result of 
the change in species 1, species 2 is now not located in a fitness peak. The landscape 
of species 2 has been modified by the adaptive motion of species 1 . 

The NKC model was analyzed computationally by Kauffman and Johnsen [37], 
and they showed that this system was dynamically very rich. They identified a 
chaotic phase, where the ecosystem is always changing and species never end up 
in a particular configuration (i.e., species do not stop at a given local peak). As we 
will discuss later, many other different models suggest that chaos can be found in 
Red Queen dynamics. They also identified a frozen phase, where all species settle 
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Fig. 5 Coevolutionary dynamics in the NKC model. Here only two species (Si and S2) are consid- 
ered, each one described by N = 3 traits, which are represented in cubic sequences spaces. Black 
circles in the hypercube nodes indicate the fitness associated to each string, and current states are 
highlighted with open circles in (a). Each bit is assumed to interact with a number of different bits 
from the other string (genome) as indicated in the lower diagrams. Once the first species changes 
by climbing to a local optimum, the landscapes of both species get modified. The second species 
now will be forced to change too, since it is now placed in a low-fit state and will next shift to 
another local peak (here indicated with a dashed circle in (b)). If no such movement is possible, 
extinction can take place. 



down to local peaks. Interestingly, for finite systems at the boundary between the 
chaotic and the frozen phase in the parameters space, small perturbations generate 
a coevolutionary avalanche of changes through the system. Typically occurring at 
critical states, the distribution of such avalanches was shown to follow a power-law. 
Kauffman and Johnsen mapped these avalanches into extinction events, suggesting 
that the number of changes in species was proportional to the extinction of less-fit 
variants. Such a result did not fit the predictions of fossil record extinctions. How- 
ever, a variation of this model by Newman and Palmer [50], which allowed changes 
in the parameters, gave an exponent which agreed with fossil record data (see also 
Section 6). 

The two phases of the NKC model can be derived from simple theoretical argu- 
ments setting K = N — 1 1 3 1 . It is known that a given species, in order to reach a local 
fitness, needs a number of walks L w , which is on average L w s» ln(N). If we assume 
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that all species are at local peaks and one of them, named a, is perturbed (i.e., is 
randomly positioned on the fitness landscape), then a will start again to climb some 
other local peak. If C is large enough (i.e., interactions among species are impor- 
tant), the other species except a can see their landscape modified also starting to 
change. Following this idea, where each adaptive walk involves a change in a given 
trait, which can in turn affect other species traits, we can determine a critical con- 
dition given by a combination of N and C able to trigger a chain reaction able to 
percolate through the system. More concretely, the probability that a given trait in a 
random species depends on species a is C/N. The critical condition is that at least 
one change in a species occurs. This actually means: 



being C c = N/ln(N). That is, when, on average, one out of C randomly chased 
genes is among the L w changed genes. In other words, when the number of traits 
is such that C > C c , interactions among different genotypes constantly modify the 
underlying fitness landscapes, scenario under which coevolutionary avalanches take 
place. 



3.4 Red queen dynamics in RNA virus 

Let us first consider a very simple example of coupled fitness landscapes and two 
identical populations climbing and competing on them. The Red Queen hypothesis 
of evolution has been widely discussed within the context of RNA viruses Il55ll70ll . 
where the dynamics of viral populations can be interpreted as a dynamical process 
of growth, competition and selection taking place in the sequence space. The fitness 
landscape for a virus is usually defined in terms of replication rate or infectivity 
or transmission. The landscape appears as a multipeaked surface, where the local 
maxima represent optimal fitness values which can be reached by mutation. Here, 
the initial condition plays an important role since depending on where the quasis- 
pecie^jis located in the sequences space, the population will evolve by exploring 
near genotypes by mutation. Competition experiments between several clonal viral 
populations Il8l l55ll provided a good illustration of two basic principles of evolution- 
ary ecology: the Red Queen dynamics and the principle of competitive exclusion. 

Experimental results were carried out with two clones of Vesicular stomatitis 
virus (VSV, see figure |6j. Such experiments involved the mixing of two clonal pop- 
ulations of VSV of equal fitness. Passage experiments allowing these populations to 
grow and compete were performed using standard virus plaque assays. More specif- 
ically, genetically marked monoclonal antibody-resistant (MARM) clones of equal 



2 The term quasispecies is used to define the heterogeneous population of viral genomes in RNA 
viruses. 
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days 

Fig. 6 Left: Vesicular stomatitis virus (VSV) virion which contains a negative-sense, single- 
stranded RNA genome. The bullet shape is characteristic of the Rhabdoviridae family (drawing 
by Ricard V. Sole). Right: Time changes of MARM-C:wt ratio in independent replicas, (a) Shows 
eighteen replicas started after passage 12. (b) Displays six replicas started from passage 14. In 
both panels the populations start to diverge after approximately 23 passages (indicated with small 
arrows). Plots obtained from |55|. 

fitness to the wild-type VSV were used and their relative frequencies were moni- 
tored along passages. 

The MARM clones only differed in a single mutation with neutral affects not 
changing viral fitness. These experiments revealed that both competing populations 
grew up showing steady increases of fitness, but, at some point, one of the two pop- 
ulations suddenly excluded the other one. The winner of this competition process 
was not always the same (see figure IS). Although the time scale of the divergence 
seemed highly predictable. The simultaneous increase of fitness of the two popu- 
lations and their predictable divergence was suggested to be a product of the Red 
Queen effect ||70ll . In the context of RNA viruses, newly arising mutants with higher 
fitness were able to outcompete lower-fitness ones. At the level of viral genomes or 
sequences, a favorable mutation within one quasispecies triggers evolutionary re- 
sponses in the second one, forcing it to evolve. Overlapped with this evolutionary 
process, and related to the dynamics, the principle of competitive exclusion is also 
at play. This principle states that when two species are strongly competing for the 
same finite resources, the fitter one asymptotically outcompetes the least fit. 

The previous experiments were modeled by Sole and collaborators |70| using 
different approaches. The simplest one was a bit-string model that considered a pop- 
ulation of N bit-strings, named S(, with sequences: Si = SjSj...Sj; with i = 1,2, ...N 
and 5/ G {0, 1}. At each time generation (passage), the algorithm repeated N times 
the following rules: a random string, say Si, of the population was chosen for repli- 
cation. Replication, proportional to replication probability r(5,-), took place by re- 
placing one of the strings of the population (also randomly chosen), say Si by a 
copy of S{. Replication presented error, at a rate (per bit and replication cycle) jj,. 
So the probability to copy exactly the same bit was 1 — ji. The mapping between 
sequence composition (genotype) and replication rate (phenotype) was done using 
the Fujiyama fitness landscape (see figure |4|b)), involving the linear relation: 
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Fig. 7 (Upper) Dynamics of the bit strings model using a population with N = 3000 strings of 
length v = 16 using fl = 10~ 3 . The main plot shows the time dynamics of both populations. The 
inset shows the mean fitness, (r), also for both populations along time (see |70|). The observed 
changes can be easily interpreted in terms of a parallel climbing of both species on their Fujiyama 
landscapes together with ongoing competition for resources. Below we illustrate this by means 
of a small, three-bit landscape. Initially both species (their populations are indicated with open 
circles) grow slowly and competition is weak. As they climb up and increase their replication 
rates, competition become strong and symmetry breaking occurs (see text). 
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As we previously explained, this fitness landscape ignores epistatic interactions. 

Simulations revealed the same behavior obtained with the experiments with VS V. 
Figure|7]displays the outcome of the model for a population of N = 3000 strings of 
length v = 16. The populations were initialized in such a way that the initial fitness 
of both populations was low, also keeping equal their mean replication rates. During 
the simulations, the strings were competing for the available space (i.e., N available 
sites). The upper panel of figure[7]shows the total population size of each population, 
which was maintained roughly constant along time. However, after approximately 
t « 20 passages, one of the two populations started outcompeting the other one, that 
finally disappeared. Once ecological competition became tight, selection pressure 
became stronger and the initial parallel growth in fitness for both populations was 
no longer sustainable. This dynamical divergence was a direct outcome of a "sym- 
metry breaking" phenomenon which explained the VSV experiments (see 0511701 
for details). 



4 Gene-for-gene and Matching Alleles Models of Coevolution 



Coevolution is the change of a biological object triggered by the change of a related 
object. Coevolution can occur at many biological scales: at the molecular level as 
correlated mutations between amino acids in a protein [92 1, or at the macroscopic 
scale as covarying traits between different interacting species in an environment. 
In coevolution, each entity exerts selective pressures on the other, thereby affecting 
each other's evolution. This process is schematically illustrated in figure[8jd). Here 
we show two fitness landscapes for preys and predators. Imagine preys are viruses 
(or cells infected by viruses) and predators are cytotoxic lymphocytes (cells of the 
immune system that kill infected cells) that act upon the activation of the adaptive 
immune response. If the virus, located in the peak with the green dot is able to mu- 
tate, visiting the highest peak, the immune system will not be able to recognize and 
remove it. However, if the virus moves towards the lower peak, which is recognized 
by dendritic cells or macrophages, able to trigger the immune response, virus pop- 
ulations with this genotype will be impaired in terms of number of particles due 
to the action of cytotoxic lymphocytes, that will remove infected cells. This simple 
example illustrates how the evolution of one of the partners influences the evolu- 
tion of the other, and viceversa, in a coevolutionary arms race. Broadly speaking, 
coevolutionary interactions can be antagonistic or mutualistic. The former involve 
negative interactions such as predation or parasitism. The latter being found when 
two or more species coevolve by means of cooperation. Coevolution can occur for 
two interacting species (pairwise coevolution) or can involve a number of different 
species, which are evolving in responses to another set of species (diffuse coevolu- 
tion). 
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Fig. 8 Coevolutionary phenomena can be interpreted as the coupling of fitness landscapes by 
means of ecological interactions (dashed arrows). We show neutral sequences spaces for predator- 
prey (host-parasite) species with perfect matching alleles interactions for (a) v = 1 , (b) v = 2 and 
(c) v = 3, being v the length of the sequences, (d) Prey's evolutionary fate will not only depend on 
its own and independent exploration of the fitness landscape, but also on predator's evolution (in- 
dicated by the small arrows). If a given prey (green circle) moves towards the highest peak, it will 
escape predator's action increasing its fitness. However, if prey climbs to the lowest peak, and the 
predator mutates moving to the same peak, host's fitness and reproductive success will diminish. 



There are at least about six proposed forms of coevolution between species, some 
involving reciprocal adaptation and others a combination of adaptation and specia- 
tion lIBTl [82l [83l . In the context of coevolution between hosts and parasites, some 
of them had a particular importance in this type of interactions. One is Ehrlich and 
Ravens [ 16] hypothesis of how the evolution of defence and counterdefence in host- 
parasite interactions may lead to the radiation of species through the process of 
escape-and-radiation coevolution. The genetic basis of infection in real ecosystems 
has been also represented by two major models: the so-called gene-for-gene (here- 
after GFG) and matching alleles (hereafter MA) models. The GFG model is based 
on data from plant-pathogen interactions, especially in the field of crop plants EOll . 
Interestingly, the first mathematical model of coevolution was explicitly based on 
assumptions of a GFG interaction |47|. Later, a multitude of real examples on GFG 
coevolution were identified between plants and pathogens, mainly between plants 
and fungi, bacteria and viruses (see [84 1 for a review). The GFG hypothesis states 
that "for each gene that conditions reaction in the host there is a corresponding gene 
that conditions pathogenicity in the parasite" ll84l[39l . The key feature of this model 
is that one parasite genotype can infect all host genotypes. 

As a difference, in the MA model, favored by invertebrate zoologists ll23l . an 
exact genetic match is required for infection (figure [8|a-c)). MA models underlie 
most of the theory constructed to understand the effects of host-parasite coevolution 
on sex and recombination ll27l[30l . Parker [52 1 pointed out the importance of MA 
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models for the study of sexual reproduction, suggesting it may hamper the generality 
of the Red Queen theory for sex. It has been argued that both GFG and MA models 
are not totally disconnected: they are two ends of a continuum (see [T| for further 
details). 

In the following sections, we will introduce and review recent models and re- 
sults about coevolving replicators with antagonistic interactions, mainly focusing in 
MA models. These models, although being suitable to analyze coevolution at small 
scales (e.g., immune system-viruses) provide good intuitions in larger scales, such 
as in spatially-extended ecosystems (Section 5.1.). In section 5.2. we will develop 
some theory aimed to describe the dynamics arising from MA predator-prey inter- 
actions. Finally, we will explore large scale coevolution by means of a complex 
network model reproducing the extinction pattern found in the fossil record data 
discussed in the Introduction Section. 



5 Minimal Coevolutionary Systems 

Coevolutionary phenomena can be studied considering minimal models. Such mod- 
els can help us to understand fundamental phenomena arising from species interac- 
tions and evolution. We notice that coevolution is a highly nonlinear phenomena, 
since species interactions give place to nonlinear couplings that can result in very 
rich and complex dynamics. In this section we will first introduce a minimal system 
of replicators with matching alleles (MA) dynamics moving, replicating and evolv- 
ing on a surface. Then, we will develop a general mathematical model describing 
MA interactions for haploid genotypes, assuming well-mixed populations thus ig- 
noring spatial correlations. As the reader will see, the dynamics of such small and 
simple systems can indeed be very complex. 



5.1 Spatial Red Queen dynamics 

At the beginning of this chapter we have illustrated the idea of coevolution with 
a very simple model simulating replicator spatial dynamics with MA interactions. 
Together with such model, other approaches have focused on the same subject by 
considering further complexity, such as spatial diffusion of replicators or larger se- 
quences spaces. Recently, Sardanyes and Sole ll59ll explored a similar system sim- 
ulating coevolution for host-parasite (prey-predator) replicators also using cellular 
automata (CA) models. 

The authors explored the spatio-temporal dynamics for three different host- 
parasite systems considering 1-bit, 2-bits and 3-bits strings [the corresponding cou- 
pled hypercubes are displayed in figure [HJa-c)]. Thus, one of the aims of their work 
was to analyze the effects of increasing the size of the sequences space in the spatio- 
temporal dynamics for MA interactions. 
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Fig. 9 Population equilibria for hosts (plain surfaces) and parasites (gridded surfaces) replicators 
in the parameters space (r^,^) for the spatially-extended model of coevolution with diffusion 
using Ji/, = [lp = 10~ 2 . Three different systems were analyzed, with: v = 1 (a), v = 2 (b) and 
v = 3 (c) (see (59j| for further details). 

The model considered two populations of replicators given by host (prey) strings 
of size v: S' h — (si ,...,s 1 ?); and by parasite (predator) strings of the same size, 
S), = (s^,...,sf) with S%,S% e {0,1} where i = l,...,N, being N the number of 
different genotypes {N = 2 V ). Both populations reproduced and evolved on a two- 
dimensional space with toroidal boundary conditions. For this model we used the 
so-called von Neumann neighborhood, which considers interactions with the four 
nearest neighbors. Specifically, the state-transition rules of the CA considered self- 
replication with mutation, decay and spatial motion of strings. Rules implemented 
are shown in Box 2. 

Box 2. State transition rules of the host-parasite CA model with matching 
allele interactions [59]: 

1 . Self-replication: If a host and a parasite occupy the same spatial position 
and have the same sequence of bits (i.e., perfect MA), the parasite elim- 
inates the host and replicates, with probability r p , to a random neighbor 
provided it is empty. If only the host lives in the cell, it replicates with 
probability r/, to a neighbor cell provided it is not occupied by another host 
string. Replication involves point mutations for host and parasites, with 
mutation probabilities /I/, and jj, p , respectively. These rules can be repre- 
sented by the following set of reactions: 

S< + tf r * (1 ~ M * )V > 25^, (3) 

S ' h + ^^%S i h +Sf. (4) 

Reactions <(3j and Q denote, respectively, error-free and erroneous host 
replication. 

&;r„(l-u„) v 

S' h +Sj, + $ M,J > 2S£, (5) 

S l h + SJ p + ^^Asj,+S l fJ. (6) 
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Similarly, reactions Q and Q represent, respectively, error-free and 
erroneous parasites replication, which is nonlinear due to the density- 
dependence of the antagonistic interaction. The parameter <5, 7 - is again the 
Kronecker 8 function where 8jj = 1 if i = j and otherwise; and # indi- 
cate some available building blocks (i.e. mononucleotides) needed to built 
new strings. The terms W-j, with k — {h, p}, correspond to the probabilities 
of erroneous replication for hosts (h) and parasites (/?), and are given by: 

being the Hamming distance between two sequences: 

^[Sf,S?] = £|4-4|. (7) 
1=1 

Equation (J7j is a function returning the number of different bits when 
comparing two sequences. Such a function is also used to determine 
the matching allele interaction between both host (S 1 -) and parasite (Sf) 
sequences, now with dH[Sf,Sf] — Yj=i \ s i ~ s f l> where s' 1 and sf represent 
the bit value (0 or 1) in the /th position in both strings. A perfect matching 
allele interaction will occur when d H [Sf,Sf] = 0. 

2. Molecular decay: Host and parasite strings decay with probability 5/, and 
8 P , respectively, according to: 

si A *, 

3. Local diffusion: Host and parasite strings move, independently and ran- 
domly, to empty neighbor cells with diffusion probabilities £)/, and D p , 
respectively. 

To simplify the analysis, all the simulations were run with maximum dif- 
fusion constants £)/, = D p — 1, also setting 5/, = 8 p = 10~ 2 . The lattice was 
randomly inoculated by either host and parasites random sequences. 



As we previously mentioned, the rules were implemented for three different sys- 
tems with different strings' lengths: v = 1, V = 2, v = 3. For all three different 
values of v, the system underwent the same three types of asymptotic dynamics: 
(i) stable coexistence of hosts and parasites with sustained fluctuations; (ii) hosts 
survival and parasites extinction; and (iii) both hosts and parasites extinction (i.e., 
coextinction). The visualization of the populations trajectories in phase space re- 
vealed the presence of chaotic coevolutionary dynamics l59l (see also l60l ). In or- 
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der to characterize the importance of these three possible asymptotic states listed 
above, we built parameter spaces considering two key evolutionary parameters of 
hosts and parasites coevolution: self-replication and mutation rates. Figure ^illus- 
trates the outcome of some simulations in the parameters space (rf,,r p ) showing the 
average population numbers for both global populations^! after some transient was 
removed for v = 1 (figure|9|a)); v = 2 (figure |9jb)); and v = 3 (figure |9jc)). A key 
result was that the increase of the length of the replicators (i.e., increase of the size 
of the sequences space) promoted stable coevolution, as shown by the reduction of 
host-parasite coextinctions (figure[9]l. 

Moreover, for each of the genotype lengths we simulated three different sce- 
narios, characterized by different values of hosts and parasites mutation rates. The 
simulations revealed that asymmetries in mutation rates between hosts and parasites 
had an important effect in the population dynamics: hosts were only able to escape 
from parasites (causing parasites extinction) if they mutated much faster. Under this 
condition, the scenario of host's survival and parasites extinction was found for ex- 
tremely low values of hosts' self -replication rates. On the contrary, when jj, p > jih, 
the region of parameters space with host and parasite extinction increased for the 
three hypercubes analyzed, indicating that when parasites evolved faster than hosts 
they were more efficient in catching hosts thus increasing coextinction phenomena. 



5.2 Dynamics of small replicators with matching-allele interactions 

The previous computational models considered antagonistic populations of bit- 
strings replicating and mutating on a surface. The same system can be investigated 
with a mathematical model by assuming no spatial correlations (i.e., infinite diffu- 
sion). A general model describing predator-prey matching-alleles (MA) interactions 
can be formulated using a time-continuous deterministic model. Assuming a perfect 
MA [see figure[8]and figure 10 a)], where each predator genotype can predate only 



on its homologous prey genotype (i.e., predator genotype i predates on prey geno- 
type i, with ; e {0, 1}), the model is given by the following system: 



X 



»•)+ 



**/, (8) 




y i = ^ Xuyi ) + tl[ £ V/ _ V( |- f /'v,, (0) 



with: 



3 by global populations we mean the sum of all possible genotypes for a given population i.e., hosts 
or parasites. 
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The state variables xi and indicate, respectively, the concentration or population 
numbers of the /th prey and of the ith predator genotype (with i = 1...2 V ), which 
define a v-dimensional sequences space J4? v . 

Note that prey genotypes have a logistic-like growth constraint in self-replication 
indicated in the first term in parenthesis of Eq. (jHJ, where £ jeJt v x j> is me total P rev 
population and J€ the carrying capacity of the prey genotypes. The logistic term 
involves exponential growth for small population numbers and saturation as popu- 
lation values approach the carrying capacity. Moreover 1 /A, is the yield coefficient 
of prey genotype i to predator genotype i. Equation (jTOJ is a Holling "type II" func- 
tional response Il6l l28ll . where predation rate is a saturating function of prey density. 
Cj and k p are constants parametrizing the saturating functional response. The con- 
stant k p , which describes the maximal predation rate, can also be interpreted as the 
predator's maximal self-replication or intrinsic growth rate. 

Both terms Y,<j> i x j ~ x u and £</>.y; —yu denote genetic diffusion by mutation 
among neighboring genotypes for both prey and predator genotypes, which are pro- 
portional to fi h and fj, p , respectively. Moreover, k 1 } denotes self-replication constant 
(intrinsic growth rates) for prey genotypes; ef and ef are decay rates which can 
be interpreted as spontaneous hydrolysis rates as well as density-independent death 
rates. If only a single genotype is present in each species, equations |8]l and Q are 
close to the well-known Rosenzweig-Mac Arthur model ll89l . 

As the reader will see, since the dimension of the dynamical system described 
by Eqs. |8]l-(|9} depends upon the length of the sequences, one may expect different 
types of dynamics for different values of v. As we will discuss in the following 
sections, were we review results for v = 1 |6D (figure 10 a)) and v = 3 l62l (figure 
[HJc)), this is the case. The fitness landscape for this predator-prey system using 
V = 1 is shown in figure [TOf a). Actually, this system was studied for two different 
fitness landscapes, given by a flat or neutral fitness landscape, with k\ = k h and 
kf = k p , V; (i.e., all genotypes share the same fitness values) and for an asymmetric 
fitness landscape with k p < K p or k p > K p (i.e., one of predator's genotypes has a 
higher fitness in terms of populations growth). For this latter case, a predator with 
a higher fitness actually means that a given genotype is more efficient in catching 
its preferred prey. For both scenarios we also assumed that mutation rates for both 
prey and predator genotypes were equal i.e. /x* = fJ. h and fif = ji p , V/ (see 16T1 for 
further details). 

The analysis of the qualitative behavior of Eqs. <[8]>-(|9j was performed assuming 
a neutral fitness landscape with ef = ef = £. This system was shown to have three 
fixed points, given by: (x* = 0,y* — 0), (xq,x\,0,0), and (xg,x\,yQ,y\). The first 
fixed point, if stable, involved predator-prey extinction. The second equilibrium, in- 
volved prey survival and predator's extinction whereas the third fixed point involved 
predator-prey coexistence, and thus it was the potential state of Red Queen dynam- 
ics. This third fixed point, under symmetry conditions, named (x*,y*), was given 
by: 
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Fig. 10 (a) Minimal predator-prey system with matching alleles interaction (vertical dashed lines) 
modeled with Eqs. {8}-{5} for v = 1. Upper and lower marbles correspond, respectively, to predator 
(y) and prey (x) genotypes, which reproduce at rates k p and k/, (circular arrows), and mutate at rates 
)lp and /!/, (straight arrows), respectively. We show an asymmetric fitness landscape for predators, 
with larger replication rates for predator genotype 1 (i.e., K p > k p ). (b) Bifurcation diagram for 
predators with genotype 0, using mutation rates (p. h and \l p are represented with thick and thin 
lines, respectively) as control parameters for the neutral fitness landscape (i.e., K p = k p ).. The 
system undergoes a Hopf bifurcation resulting in a permanent oscillatory behavior governed by a 
periodic orbit 1 6 1 1 . 
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To study the stability of this fixed point it was further assumed that k h = k p = 1, 
also considering symmetry in decay rates i.e., e h = £ p = £. Under this conditions, 
the fixed point reads: 



x = 



y 



l-e' 

(l-4£ + £ 2 ) 



(l-e) 2 ' 

After some algebra, and after fixing jj, h = jj. p = fi, a critical mutation rate causing a 
Hopf bifurcation was identified at: 

_ e(l -4e + £ 2 ) 
4(1 -e) • 

Such a bifurcation, which involves the creation of a periodic orbit causing sus- 



tained, periodic oscillations, was confirmed by numerical simulations (figure 10 b)). 
Interestingly, the same bifurcation was also numerically found for the fitness land- 
scape with asymmetries in predator's replication rates. Counterintuitively, the asym- 
metric fitness landscape revealed that the most efficient predator genotypes achieved 
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lower population equilibria (see [61 1 for details). Our results identifying periodic 
Red Queen dynamics were in agreement with other mathematical models on revo- 
lution (see [11]) . 



5.3 Chaotic Red Queen attractors 

So far, we have discussed the dynamical behavior of Eq. |8]l-(|9) for v = 1, which can 
be governed by stable fixed points as well as by a periodic orbit causing sustained 
and regular oscillations of predator-prey genotypes populations. The same model, 
analyzed for v = 3, revealed much richer dynamics: under some parameter regions, 
both populations behave chaotically. Hence, similarly to what is known as diffusion- 
induced chaos [53 1, it was found that the simplest system (with V = 1), governed by 
a periodic orbit, could be governed by chaotic attractors at increasing the number of 
available alleles (more available nodes in sequences space). 

It is known that dynamical systems governed by a periodic orbit can become un- 
stabilized to chaos when spatial correlations and diffusion are included. Actually, 
some decades ago a great deal of attention was paid to self-organization processes 
in reaction-diffusion systems, and their relevance in chemistry, physics and biology 
was repeatedly stressed [51, 24, 88]. In this sense, numerical investigations of the 
spatially-extended Belousov-Zhabotinsky chemical reaction showed the presence of 
chaotically oscillating structures. Moreover, diffusion-induced chaos has also been 
discussed in the context of spatial ecological dynamics ll53l . We notice that MA 
chaotic dynamics can be also interpreted from the perspective of dynamical unsta- 
bilization due to diffusion in space. That is, each sequence of the sequences space 
can be interpreted as a patch, and populations can diffuse between patches because 
of mutation (i.e., diffusion in sequences space). Under this view, the oscillatory be- 
havior of variables Xj and y, in Eqs. ([8]l-(|9]) for v = 1, becomes unstabilized to chaos 
for v = 3. 



Figure 1 1 shows the chaotic coevolutionary dynamics for the 3-bits sequences 
modeled with Eqs. ([8])-(|9]l. The time series in figure [TTJa) are represented for the 
global populations of predators (y) and preys (x). In figure [TTJ a) and (b) we show, 
respectively, the chaotic attractors for global host-parasite populations (represented 
in the phase space (x,y)) as well as the attractor for parasite genotypes 010, 111, 
and 100 (see ll62ll for further details). 

Our previous results suggested that Red Queen dynamics can be chaotic even for 
small haploid replicator systems with MA interactions. It was previously shown that 
large networks with host-parasitoid replicators can also behave chaotically. More 
specifically, Kaneko and Ikegami ||3Tll36l characterized the so-called homeochaos in 
multi-species models with antagonistic interactions and evolution. They suggested 
that chaos, more than a destabilizing behavior 0, could involve stability in multi- 
species ecosystems through a weak, chaotic state arising in high-dimensional dy- 
namical systems. Roughly, homeochaos was suggested to suppress strong chaos 
causing large fluctuations that could near populations to extinction. Homeochaos is 
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Fig. 11 Red Queen chaotic dynamics for Eqs. |8)-((9j using v = 3. (a) Global population time dy- 
namics of parasites (y) and host (x). (b) Strange attractor governing host-parasite global dynamics. 
In (c) we show the chaotic attractor in the parasites three-dimensional phase space for genotypes 
(OlOy, lily, lOOy) (see l62l l. 

characterized by many positive, but close to zero Lyapunov exponents (i.e., a type 
of hyperchaos). Such a property of the spectrum of Lyapunov exponents involves 
narrow chaotic fluctuations with small amplitude, which are able to keep population 
numbers far away from attractors involving extinction. The concept of homeochaos 
was later extended to low-dimensional systems, and its role was discussed in both 
deterministic and stochastic host-parasitoid models with discrete time generations 

m. 

Chaotic evolutionary dynamics have been found in other theoretical studies 
of genetic polymorphisms under frequency-dependent selection (see for example 
B31l67l[T9 l). Moreover, Dercole and colleagues 1111 recently showed that predator- 
prey coevolutionary models governed by periodic fluctuations became chaotic when 
the system is embedded in a three-species food chain model by the addition of a 
superpredator able to coevolve. These authors argued that over space, genetically- 
driven chaos may cause evolutionary divergence of local metapopulations, even un- 
der the absence of environmental change, thus promoting genetic diversity among 
ecological communities over long evolutionary time. 
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6 Large Scale Coevolution on Complex Networks 



Our last example deals with the large-scale evolutionary dynamics. The study of 
large scale coevolution was also performed in multi-species models using complex 
networks theory. A very simple model of large-scale evolution involving a set of 
N interacting species can be easily defined l72l l73l l74l l75l . In this model, species 
interactions are introduced by means of a N x N connectivity matrix W = (Wij). 
Evolution for this system is introduced through changes in its elements. Similarly 
to some of the models previously presented, here the "state" of each species is de- 
scribed by a binary variable 5, E {0, 1} (z —1,2, ...,N), for the z'-th species, with 
5, = 1 or Sj = if the species is alive or extinct, respectively. So the whole ecosys- 
tem is described in terms of a simple directed graph where the connections are ini- 
tially set to random values. Each species receives either positive or negative inputs 
from other species. These signs indicate that the given species is favored or harmed 
by the species which send the input. For instance, a negative input would correspond 
to the interaction with a predator or with a parasite. Alternatively, a positive input 
would correspond to mutualism or symbiosis. Such a model, in its simplest form, 
can be formulated in terms of a set of rules displayed in Box 3. 



Box 3. Rules of the network coevolution model (illustrated in figure 12 a-d)): 



1. Random changes in the connectivity matrix. At each generation, we 
select one input connection for each species and assign it to a new, 
random value without regard for the previous state of the connection. This 
rule introduces changes into the web, which can be due to evolutionary 
responses or to environmental changes of some sort. In other words, 
changes derived from coevolution among two species, innovation at the 
species level and/or environmental-driven changes are lumped together 
within this rule. 



Extinction. Changes in the connectivity will eventually lead to extinctions. 
Extinction events are decided by computing the total sum of the inputs 
for each species. This sum, if negative will involve the extinction of the 
species (Si — 0), and all its connections are removed. Otherwise, nothing 
happens (5, = 1). Hence, the state of the z-th species is updated following 
the following dynamical equation: 



i>(z,j) 



Sj(f + 1) = * 
where <P(z) = 1 for positive z and zero otherwise. 



3. Diversification. A number of species can disappear due to the extinction 
rule, leaving empty sites. These sites will be refilled by diversification: each 
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extinct species is replaced by a randomly chosen survivor. The replacement 
is made by simply copying the connections of the survivor into the empty 
site. 



This model shows a strongly nonlinear behavior with avalanches of extinction 
as well as the correct power law distribution of extinction sizes l72l l74l l44l |78l 
(although with an exponent typically close but higher than a — 2, see also |15|). 
The outcome of the model was that both small and very large events were generated 
by the same dynamical rules. Most of the times, the extinction of a given species 
had no consequences for the other species. But from time to time, a given (keystone) 
species with positive inputs to others disappeared. The removal of this species was 
suggested to have a destabilizing effect on others, able to cause further propagation 
leading to mass extinction events. 

What is the origin of such extinction patterns? We first need to see how a given 
species can shift towards a negative sum of inputs. The reason is easily understand- 
able from rule 1 above. 

Since changes of links among species are random and the new values are chosen 
from a uniform distribution an expected consequence is that, in the long run, the 
sum of inputs will decay to zero. If we look at the sign of the links, so that the 
probability of finding positive links, P(W + ) = P[Wjj > 0]; and the probability of 
finding negative links P(W~) — P[Wij < 0]. The time evolution of the positive links 
can be described in terms of a master equation, given by: 

dP(W + t) 

— v — = w(W~ -> W + )P(W+) - w(W + -> W~)P(W~), 

at 

where P(W + ) +P(W~) = 1, starting for example from an initial condition P(W + ,0) - 
P . Since, from the first rule, we have w(W~ -t W + ) = w(W + -> W~) = 1 /IN, the 
master equation reads: 

and an exponential decay is obtained: 



P(W + ,t)= 1 -[l + (2P Q -l)e-'/ N 



As a consequence, the sum of inputs J?, = E/Wy will also decay exponentially: 
&i(t) ~ e~'/ N , predicting an exponential decay in the probabilities of survival, as 
expected from the Red Queen hypothesis (see Introduction). This rule actually in- 
troduces the basic ingredients of Van Valen's theory. All species in the system keep 
changing all the time (either due to biological or environmental causes) eventually 
reaching extinction. The ultimate fate of all species is to get extinct, and so an ex- 
ponential decay in the survival probability will be observed. Here, however, there is 
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Fig. 12 (a-d) Rules of the evolution model exemplified with a network with N = 7 nodes and a 
given initial connectivity (a). Negative and positive interactions are indicated with black and white 
arrows, respectively. The first step of the model is to modify the links (b): here different pairs 
of interactions are found, including mutualism, parasitism, predation and competition. For this 
particular case, species 3 and 4 become nonviable, dying in the next step (c). Species 6 is selected 
(among the survivors) and copied into the two empty sites (Rule 3). Such copies carry the links 
of their parent species, which will be modified again by the first rule in an iterative process. This 
set of rules generates a very complex dynamical pattern of species evolution. In (e) we display an 
example of the survival curves obtained from our model (compare with figure [TJ. In (f) we also 
show the evolution of the local fields over time. A mixture of slow and rapid changes occur, in a 
punctuated fashion (see 1 73 1 for further details). 



no intrinsic, species-level variability (in terms of a genotype) and the fate of a given 
species will be dominated by network responses and chance. Ecological-driven phe- 
nomena are the key forces in the long run, although small-scale events are taking 
place all the time. The nature of the decay turns to be exponential on average but 
episodic when looking at pseudocohorts (see figure [TJ, consistently with Raup's 
analysis commented on in the Introduction of this chapter (compare figures [T] and 

An analytic study of the previous model is rather difficult because of the random 
nature of the interaction matrix. The model can, however, be simplified by mapping 
the set of rules into a linear model ]44ll78l . giving place to a mean field approach to 
the network model (see Box 4). 

Box 4. Consider a set of N species, characterized by a single integer quantity 
0; (i = 1,2,...,N). This quantity will play the role of the internal field. Each 
species is now represented by this single (integer) number 0,- € {— N, — N + 
1,..., — 1,0, 1,...,N— l,N}, which represents the sum of inputs from other 
species. The dynamics consists of three steps: (a) with probability P — 1/2, 
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0? - > <Pi — 1> otherwise no change occurs (this is equivalent to the randomiza- 
tion rule in the network model); (b) all species with 0,- < <j> c (below a given 
threshold) are extinct. Here we use <j> c = but other choices give the same 
results. The number of extinct species, < E < N, gives the size of the ex- 
tinction event. All E extinct species are replaced by survivors. Specifically, for 
each extinct site (i. e. when 0; < C ) we choose one of the N — E survivors 0^ 
and update ; - to: ; - = 0^; (c) after an extinction event, a wide reorganization 
of the web structure occurs [75]. In this simplified model this is introduced as 
a coherent shock. Each of the survivors are updated as 0^ = 0^ + q(E), where 
q(E) is a random integer between — E and +E. This mean-field approach de- 
fines a three-step process. If N((j>) indicates the frequency of species having a 
local field 0, we have: 

jV(0,r + l/3) = ^V(0, f) + ^(0 + l, t), 

AT(0,f + 2/3)=iV(0,/ + l/3)+iV(0,/ + l/3)^^ wj P(m), 
if > and zero otherwise. Finally: 
NQ,t + 1) =N(#,t + 2/3) -NQ,t +1/3) + £ N(<l> +q,t + l/3)P(q), 

q>-(j> 

from these equations, the full master equation for the dynamics reads: 

N(<j)+q,t)-N{<t>+q+l,t) 



1 , , , , w-i mP(m) 

Where two basic statistical distributions, which are self-consistently related, 
have been used. These are: 

P*(q)=Z^e(m-\q\), 

which is an exact equation giving the probability of having a shock of size 
q. The second is P e (m), the extinction probability for an event of size m. We 
have a mean-field approximation relating both distributions: 

r ?-i 

P e (m)=Z,P*(q)S £>(0)-m 

9 >=1 
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The last equation introduces the average profile N((j>), i. e. the time-averaged 
distribution of 0-values. For the mesoscopic regime 1 ^> q ^> N, by applying 
a Taylor expansion to the master equation, i. e. 



1 +^ P( m ) f A . . dN 




and using a continuous approximation, it is easy to see that the previous equa- 
tion reads: 



dm 



P(m) \ [ N(<j)+q) 
2m | L N(<j>) 



dLnN 



d(j) 



1 dLnN 
2 \ 2+ ^V 



mP(m) 
N — m' 



dm = 0. 



Assuming that N((j>) decays exponentially, i. e. N(<p) = exp(—c<j)/N), we 
can integrate each part of the last equation, using N(<j> + q)/N((j>) = 
exp(—cq/N) rj 1 —cq/N. The first term cancels exactly, the second gives 
—2c/N and the third scales as (1 — 0(1 /N))N l ~ T . So the previous equation 
leads to: 

1 

,N . 



2c 

'n" 



-N l - T G 



1-0 



= 0, 



in order to satisfy this equality, we have t = 2, which gives us the scaling 
exponent for the extinction distribution. Hence, in agreement with Burlando's 
analysis, the taxonomy that emerges from this model also displays fractal be- 
havior (with an exponent aj « 2). 



These models, able to reproduce observed patterns, can have important implica- 
tions for evolutionary theory. An intense debate over the last decades has concerned 
the basic mechanisms operating at different temporal scales. Some authors (spe- 
cially in the field of population genetics) suggested that the rules operating at small 
scales (i.e., microevolutionary events) can be directly translated into the process 
of macroevolution lf25l [T4l l42l . However others authors like Stephen Jay Gould, 
claimed that different processes are at work in evolution at different scales 11221 . al- 
though no well-defined mechanism for such decoupling was proposed. We want to 
notice that the network organization of ecologies, changing in a coevolving land- 
scape, suggests a possible source of decoupling the micro- and macro scales. More- 
over, these models can also help understanding the complex dynamical behavior of 
large extinctions and their aftermath 1176117711581 171. 
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7 Conclusions 



Coevolutionary dynamics introduces an additional layer beyond single-species evo- 
lution. Coevolution pervades biology on multiple scales but its role and impact -as 
illustrated by our previous examples- is rather different at each scale. In microor- 
ganisms, the changes that couple diverse species (as hosts and pathogens) within a 
given ecosystem are associated to many different molecular events related to mem- 
brane receptors, production and pumping of toxins, development of aggregates or 
resistance to antibiotics, to cite just a few. Arm races are known to occur and take 
place on short time scales. 

Coevolution occurs in other systems, including technological ones. Coevolution 
between predators and prey predate at least part of the evolutionary events that trig- 
gered the emergence of complex animals at the base of the Cambrian explosion. 
It is likely that the so called Ediacaran fauna, dominated by simple, filtering or- 
ganisms with small developmental complexity became replaced by the well known, 
Burguess-Shale pattern as a consequence of predator-prey arm races. Many chal- 
lenges lie ahead in our understanding of how coevolution shaped biological com- 
plexity and how to properly approach it from a theoretical perspective. Among 
other questions, we still need to understand how to connect ecological networks 
and coevolving landscapes, how to place these landscapes in the middle of the mul- 
tidimensional space involving development, ecology and the environment, and what 
universal trends are to be found in their structure and dynamics. The previous exam- 
ples only provide a glimpse of the richness and complexity, but they also illustrate 
the power of simple models able to address relevant questions. 
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